/*
** ILevinson.cs
**
**
** Copyright (c) 2003-2006, dnAnalytics. All rights reserved.
*/

#region Using directives

using System;
using dnAnalytics.Math;
using dnAnalytics.Exceptions;
using dnAnalytics.LinearAlgebra;
using dnAnalytics.Resources;

#endregion Using directives

namespace dnAnalytics.LinearAlgebra.Decomposition
{
    /// <summary>
    /// A Levinson decomposition for square Toeplitz systems.
    /// </summary>
    /// <remarks>
    /// This class provides members for inverting the Toeplitz matrix (see <see cref="Inverse()"/> member) and
    /// calculating the determinant of the matrix (see <see cref="Determinant"/> member).
    /// <para>
    /// The class implements a <b>UDL</b> decomposition of the inverse of the 
    /// square Toeplitz matrix. The decomposition is based upon Levinson's algorithm. As
    /// a consequence, all operations require approximately <b>N</b> squared FLOPS, where
    /// <b>N</b> is the matrix order. This is significantly faster than <b>LU</b> factorization,
    /// which requires <b>N</b> cubed FLOPS.
    /// </para>
    /// <para>
    /// One disadvantage of Levinson's algorithm is that all the leading sub-matrices
    /// of the Toeplitz matrix must be non-singular. During the decomposition,
    /// sub-matrices are checked to ensure that they are non-singular. When a singular
    /// matrix is found the decomposition is halted and an internal flag is set. The
    /// <see cref="IsSingular"/> property may be used to access the flag, to determine
    /// if any singular matrices were detected.
    /// </para>
    /// <para>
    /// A outline of this approach to the UDL decomposition of inverse Toeplitz matrices is
    /// found in the following reference:
    /// </para>
    /// <para>
    /// <i>Sun-Yuan Kung and Yu Hen Hu</i>, A Highly Concurrent Algorithm and Pipelined
    /// Architecture for Solving Toeplitz Systems, IEEE Transactions on Acoustics,
    /// Speech and Signal Processing, Volume ASSP-31, Number 1, February 1983, pages
    /// 66 - 75.
    /// </para>
    /// <para>
    /// This class is optimised for non-symmetric Toeplitz matrices. If you wish
    /// to solve a symmetric Toeplitz system, use the <see cref="SymmetricLevinson"/> class,
    /// which is twice as fast as this class.
    /// </para>
    /// </remarks>
    public sealed class Levinson : Algorithm
    {
        #region Fields
        /// <summary>
        /// The left-most column of the Toeplitz matrix.
        /// </summary>
        private readonly double[] leftColumn;

        /// <summary>
        /// The top-most row of the Toeplitz matrix.
        /// </summary>
        private readonly double[] topRow;

        /// <summary>
        /// The order of the square Toeplitz matrix
        /// </summary>
        private readonly int order;

        /// <summary>
        /// Flag to indicate if the matrix or one of the sub-matrices is singular.
        /// </summary>
        private bool isSingular;

        /// <summary>
        /// The lower left triangular matrix of the factorisation.
        /// </summary>
        private double[][] lowerTriangle;

        /// <summary>
        /// The leading diagonal of the diagonal matrix of the factorisation.
        /// </summary>
        private double[] diagonal;

        /// <summary>
        /// The upper right triangular matrix of the factorisation.
        /// </summary>
        private double[][] upperTriangle;

        private double det = double.MinValue;

        #endregion Fields

        #region Property
        /// <summary>
        /// Get the order of the Toeplitz matrix.
        /// </summary>
        public int Order
        {
            get
            {
                return order;
            }
        }
        #endregion Property

        #region Constructors

        /// <summary>
        /// Constructor with <see cref="Vector"/> parameters.
        /// </summary>
        /// <param name="leftColumn">The left-most column of the Toeplitz matrix.</param>
        /// <param name="topRow">The top-most row of the Toeplitz matrix.</param>
        /// <exception cref="ArgumentNullException">
        /// <paramref name="leftColumn"/> is a null reference,
        /// <para>or</para>
        /// <para><paramref name="topRow"/> is a null reference.</para>
        /// </exception>
        /// <exception cref="RankException">
        /// The length of <paramref name="leftColumn"/> is zero,
        /// <para>or</para>
        /// <para>the length of <paramref name="leftColumn"/> does not match that of <paramref name="topRow"/>.</para>
        /// </exception>
        /// <exception cref="ArithmeticException">
        /// The values of the first element of <paramref name="leftColumn"/> and <paramref name="topRow"/> are not equal.
        /// </exception>
        public Levinson(Vector leftColumn, Vector topRow)
        {
            // check parameters
            if (leftColumn == null)
            {
                throw new System.ArgumentNullException("leftColumn", Strings.NullParameterException);
            }
            if (leftColumn.Count == 0)
            {
                throw new RankException(Strings.LengthOfZero);
            }
            if (topRow == null)
            {
                throw new System.ArgumentNullException("topRow", Strings.NullParameterException);
            }
            if (leftColumn.Count != topRow.Count)
            {
                throw new RankException(Strings.ParametersNotConformable);
            }
            if (leftColumn[0] != topRow[0])
            {
                throw new ArithmeticException(Strings.FirstElementsNotEqual);
            }
            // save the vectors
            this.leftColumn = leftColumn.ToArray();
            this.topRow = topRow.ToArray();
            order = leftColumn.Count;

            CreateArrays();
        }

        /// <summary>
        /// Constructor with <b>double[]</b> parameters.
        /// </summary>
        /// <param name="leftColumn">The left-most column of the Toeplitz matrix.</param>
        /// <param name="topRow">The top-most row of the Toeplitz matrix.</param>
        /// <exception cref="ArgumentNullException">
        /// <paramref name="leftColumn"/> is a null reference,
        /// <para>or</para>
        /// <para><paramref name="topRow"/> is a null reference.</para>
        /// </exception>
        /// <exception cref="RankException">
        /// The length of <paramref name="leftColumn"/> is zero,
        /// <para>or</para>
        /// <para>the length of <paramref name="leftColumn"/> does not match that of <paramref name="topRow"/>.</para>
        /// </exception>
        /// <exception cref="ArithmeticException">
        /// The values of the first element of <paramref name="leftColumn"/> and <paramref name="topRow"/> are not equal.
        /// </exception>
        public Levinson(double[] leftColumn, double[] topRow)
        {
            // check parameters
            if (leftColumn == null)
            {
                throw new System.ArgumentNullException("leftColumn", Strings.NullParameterException);
            }
            else if (leftColumn.Length == 0)
            {
                throw new RankException(Strings.LengthOfZero);
            }
            else if (topRow == null)
            {
                throw new System.ArgumentNullException("topRow", Strings.NullParameterException);
            }
            else if (leftColumn.Length != topRow.Length)
            {
                throw new RankException(Strings.ParametersNotConformable);
            }
            else if (leftColumn[0] != topRow[0])
            {
                throw new ArithmeticException(Strings.FirstElementsNotEqual);
            }

            // save the vectors
            this.leftColumn = (double[])leftColumn.Clone();
            this.topRow = (double[])topRow.Clone();
            order = leftColumn.Length;

            CreateArrays();
        }
        #endregion Constructors

        #region Public Members
        /// <summary>
        /// Check if the Toeplitz matrix or any leading sub-matrices are singular.
        /// </summary>
        /// <returns><b>true</b> if the Toeplitz matrix or any leading sub-matrices are singular; otherwise <b>false</b>. </returns>
        /// <remarks>
        /// If the Toeplitz matrix or any leading sub-matrices are singular, it is
        /// not possible to complete the <b>UDL</b> decomposition using the Levinson
        /// algorithm.
        /// </remarks>
        public bool IsSingular()
        {
            Compute();
            return isSingular;
        }

        ///<summary>
        /// Get the lower triangle matrix of the UDL factorisation.
        /// </summary>
        /// <returns>The lower triangle matrix of the UDL factorisation.</returns>
        /// <remarks>
        /// It is recommended that the <see cref="IsSingular"/> property
        /// be checked to see if the decomposition was completed, before attempting
        /// to obtain the lower triangle matrix.
        /// </remarks>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        public Matrix Lower()
        {
            Compute();

            // check if there was a result
            if (isSingular)
            {
                throw new SingularMatrixException();
            }

            Matrix result = MatrixBuilder.CreateMatrix(order, MatrixType.Dense);
            Lower(result);
            return result;
        }

        ///<summary>
        /// Get the lower triangle matrix of the UDL factorisation.
        /// </summary>
        /// <param name="result">The <see cref="Matrix"/> to copy the lower triangle matrix 
        /// of the UDL factorisation into.</param>
        /// <remarks>
        /// It is recommended that the <see cref="IsSingular"/> property
        /// be checked to see if the decomposition was completed, before attempting
        /// to obtain the lower triangle matrix.
        /// </remarks>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException">If <paramref name="result"/> is not a square 
        /// <see cref="Matrix"/> with the correct order.</exception>
        public void Lower(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if (result.Rows != order || result.Columns != order)
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }

            Compute();

            // check if there was a result
            if (isSingular)
            {
                throw new SingularMatrixException();
            }

            for (int i = 0; i < order; i++)
            {
                for (int j = 0; j <= i; j++)
                {
                    result.ValueAt(i, j, lowerTriangle[i][j]);
                }
            }
        }

        ///<summary>
        /// Get the diagonal matrix of the UDL factorisation.</summary>
        /// <remarks>
        /// <returns>The diagonal matrix of the UDL factorisation.</returns>
        /// It is recommended that the <see cref="IsSingular"/> property
        /// be checked to see if the decomposition was completed, before attempting
        /// to obtain the diagonal matrix.
        /// </remarks>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        public Matrix Diagonal()
        {
            Compute();
            // check if there was a result
            if (isSingular)
            {
                throw new SingularMatrixException();
            }

            Matrix result = MatrixBuilder.CreateMatrix(order, MatrixType.Dense);
            Diagonal(result);
            return result;
        }

        ///<summary>
        /// Get the diagonal matrix of the UDL factorisation.</summary>
        /// <param name="result">The <see cref="Matrix"/> to copy the diagonal matrix 
        /// of the UDL factorisation into.</param>
        /// <remarks>
        /// It is recommended that the <see cref="IsSingular"/> property
        /// be checked to see if the decomposition was completed, before attempting
        /// to obtain the diagonal matrix.
        /// </remarks>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException">If <paramref name="result"/> is not a square 
        /// <see cref="Matrix"/> with the correct order.</exception>
        public void Diagonal(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if (result.Rows != order || result.Columns != order)
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }

            Compute();

            // check if there was a result
            if (isSingular)
            {
                throw new SingularMatrixException();
            }

            result.SetDiagonal(diagonal);
        }

        ///<summary>
        /// Get the upper triangle matrix of the UDL factorisation.</summary>
        /// <remarks>
        /// <returns>The upper triangle matrix of the UDL factorisation.</returns>
        /// It is recommended that the <see cref="IsSingular"/> property
        /// be checked to see if the decomposition was completed, before attempting
        /// to obtain the upper triangle matrix.
        /// </remarks>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        public Matrix Upper()
        {
            Compute();

            // check if there was a result
            if (isSingular)
            {
                throw new SingularMatrixException();
            }

            Matrix result = MatrixBuilder.CreateMatrix(order, MatrixType.Dense);
            Upper(result);
            return result;
        }

        /// <summary>
        /// Get the upper triangle matrix of the UDL factorisation.</summary>
        /// <param name="result">The <see cref="Matrix"/> to copy the upper triangle 
        /// matrix of the UDL factorisation into.</param>
        /// <remarks>It is recommended that the <see cref="IsSingular"/> property
        /// be checked to see if the decomposition was completed, before attempting
        /// to obtain the upper triangle matrix.
        /// </remarks>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException">If <paramref name="result"/> is not a square 
        /// <see cref="Matrix"/> with the correct order.</exception>
        public void Upper(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if (result.Rows != order || result.Columns != order)
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }

            Compute();

            // check if there was a result
            if (isSingular)
            {
                throw new SingularMatrixException();
            }

            for (int i = 0; i < order; i++)
            {
                for (int j = 0; j <= i; j++)
                {
                    result.ValueAt(j, i, upperTriangle[i][j]);
                }
            }
        }

        /// <summary>
        /// Get a copy of the Toeplitz matrix.
        /// </summary>
        /// <returns>A copy of the Toeplitz matrix.</returns>
        public Matrix ToeplitzMatrix()
        {
            Compute();

            // check if there was a result
            if (isSingular)
            {
                throw new SingularMatrixException();
            }

            Matrix result = MatrixBuilder.CreateMatrix(order, MatrixType.Dense);
            ToeplitzMatrix(result);
            return result;
        }

        /// <summary>
        /// Get a copy of the Toeplitz matrix.
        /// </summary>
        /// <param name="result">The <see cref="Matrix"/> to copy the Toeplitz matrix into.</param>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException">If <paramref name="result"/> is not a square 
        /// <see cref="Matrix"/> with the correct order.</exception>
        public void ToeplitzMatrix(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if (result.Rows != order || result.Columns != order)
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }

            Compute();

            // check if there was a result
            if (isSingular)
            {
                throw new SingularMatrixException();
            }

            int i, j, k;

            // fill lower triangle
            for (i = 0; i < order; i++)
            {
                for (j = i, k = 0; j < order; j++, k++)
                {
                    result.ValueAt(j, i, leftColumn[k]);
                }
            }

            // fill upper triangle
            for (i = 0; i < order; i++)
            {
                for (j = i + 1, k = 1; j < order; j++, k++)
                {
                    result.ValueAt(i, j, topRow[k]);
                }
            }
        }

        /// <summary>
        /// Get the determinant of the Toeplitz matrix.
        /// </summary>
        /// <returns>The determinant.</returns>
        /// <remarks>
        /// It is recommended that the <see cref="IsSingular"/> property
        /// be checked to see if the decomposition was completed, before attempting
        /// to obtain the determinant.
        /// </remarks>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        public double Determinant()
        {
            Compute();
            // check if there was a result
            if (isSingular)
            {
                return 0;
            }
            if (det == double.MinValue)
            {
                det = 1.0;
                for (int i = 0; i < order; i++)
                {
                    det *= diagonal[i];
                }
                det = 1.0 / det;
            }
            return det;
        }

        /// <summary>
        /// Get the inverse of the Toeplitz matrix.
        /// </summary>
        /// <returns>The inverse matrix.</returns>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        /// <remarks>
        /// The class implicitly decomposes the inverse Toeplitz matrix into a <b>UDL</b> factorisation
        /// using the Levinson algorithm, before using an extended version of Trench's algorithm to complete
        /// the inversion.
        /// <para>
        /// The extended version of Trench's algorithm requires approximately <b>N</b> squared FLOPS,
        /// compared to <b>N</b> cubed FLOPS if we simply multiplied the <b>UDL</b> factors
        /// (<b>N</b> is the matrix order).
        /// </para>
        /// </remarks>
        public Matrix Inverse()
        {
            Compute();

            if (isSingular)
            {
                throw new SingularMatrixException();
            }

            Matrix result = MatrixBuilder.CreateMatrix(order, MatrixType.Dense);
            Inverse(result);
            return result;

        }

        /// <summary>
        /// Get the inverse of the Toeplitz matrix.
        /// </summary>
        /// <param name="result">The <see cref="Matrix"/> to copy the inverse of the Toeplitz matrix into.</param>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        /// <exception cref="ArgumentNullException">If <paramref name="result"/> is <b>null</b>.</exception>
        /// <exception cref="NotConformableException">If <paramref name="result"/> is not a square 
        /// <see cref="Matrix"/> with the correct order.</exception>
        /// <remarks>
        /// The class implicitly decomposes the inverse Toeplitz matrix into a <b>UDL</b> factorisation
        /// using the Levinson algorithm, before using an extended version of Trench's algorithm to complete
        /// the inversion.
        /// <para>
        /// The extended version of Trench's algorithm requires approximately <b>N</b> squared FLOPS,
        /// compared to <b>N</b> cubed FLOPS if we simply multiplied the <b>UDL</b> factors
        /// (<b>N</b> is the matrix order).
        /// </para>
        /// </remarks>
        public void Inverse(Matrix result)
        {
            if (result == null)
            {
                throw new ArgumentNullException("result", Strings.NullParameterException);
            }

            if (result.Rows != order || result.Columns != order)
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }

            Compute();

            if (isSingular)
            {
                throw new SingularMatrixException();
            }

            double[] A = lowerTriangle[order - 1];
            double[] B = upperTriangle[order - 1];
            double A1, B1, scale, temp;
            int i, j, k, l;

            // setup the first row
            scale = diagonal[order - 1];
            for (i = 0, j = order - 1; i < order; i++, j--)
            {
                result.ValueAt(0, i, scale * B[j]);
            }

            // calculate values in the rest of the wedge
            for (i = 1; i < order; i++)
            {
                A1 = scale * A[order - i - 1];
                B1 = scale * B[i - 1];
                result.ValueAt(i, 0, A1);
                for (j = 1, k = 0, l = order - 2; j < order - i; j++, k++, l--)
                {
                    temp = result.ValueAt(i - 1, k) + A1 * B[l] - B1 * A[k];
                    result.ValueAt(i, j, temp);
                }
            }

            // inverse is a persymmetric matrix.
            for (i = 0, j = order - 1; i < order; i++, j--)
            {
                for (k = 0, l = order - 1; k < j; k++, l--)
                {
                    result.ValueAt(l, j, result.ValueAt(i, k));
                }
            }
        }

        /// <summary>
        /// Solve a square Toeplitz system with a right-side vector.
        /// </summary>
        /// <param name="Y">The right-hand side of the system.</param>
        /// <returns>The solution vector.</returns>
        /// <exception cref="ArgumentNullException">
        /// If <paramref name="Y"/> is <b>null</b>.
        /// </exception>
        /// <exception cref="NotConformableException">
        /// The length of <paramref name="Y"/> is not equal to the number of rows in the Toeplitz matrix.
        /// </exception>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        /// <remarks>
        /// This member solves the linear system <b>TX</b> = <b>Y</b>, where <b>T</b> is
        /// the square Toeplitz matrix, <b>X</b> is the unknown solution vector
        /// and <b>Y</b> is a known vector.
        /// <para>
        /// The class implicitly decomposes the inverse Toeplitz matrix into a <b>UDL</b> factorisation
        /// using the Levinson algorithm, before calculating the solution vector.
        /// </para>
        /// </remarks>
        public Vector Solve(Vector Y)
        {
            // check parameters
            if (Y == null)
            {
                throw new System.ArgumentNullException("Y", Strings.NullParameterException);
            }
            if (order != Y.Count)
            {
                throw new NotConformableException("Y", Strings.ParameterNotConformable);
            }

            if (IsSingular())
            {
                throw new SingularMatrixException();
            }

            Vector result = VectorBuilder.CreateVector(order, VectorType.Dense);
            Solve(Y, result);
            return result;
        }

        /// <summary>
        /// Solve a square Toeplitz system with a right-side vector.
        /// </summary>
        /// <param name="Y">The right-hand side of the system.</param>
        /// <param name="result">The solution vector.</param>
        /// <exception cref="ArgumentNullException">
        /// If either <paramref name="Y"/> or <paramref name="result"/> is <b>null</b>.
        /// </exception>
        /// <exception cref="NotConformableException">
        /// If the length of <paramref name="Y"/> or <paramref name="result"/> is not equal to the number of rows in the Toeplitz matrix.
        /// </exception>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        /// <remarks>
        /// This member solves the linear system <b>TX</b> = <b>Y</b>, where <b>T</b> is
        /// the square Toeplitz matrix, <b>X</b> is the unknown solution vector
        /// and <b>Y</b> is a known vector.
        /// <para>
        /// The class implicitly decomposes the inverse Toeplitz matrix into a <b>UDL</b> factorisation
        /// using the Levinson algorithm, before calculating the solution vector.
        /// </para>
        /// </remarks>        
        public void Solve(Vector Y, Vector result)
        {
            // check parameters
            if (Y == null)
            {
                throw new System.ArgumentNullException("Y", Strings.NullParameterException);
            }
            if (order != Y.Count)
            {
                throw new NotConformableException("Y", Strings.ParameterNotConformable);
            }

            // check parameters
            if (result == null)
            {
                throw new System.ArgumentNullException("result", Strings.NullParameterException);
            }
            if (order != result.Count)
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }

            if (IsSingular())
            {
                throw new SingularMatrixException();
            }

            // allocate memory for solution
            double inner;
            int i, j;
            double[] A, B;

            for (i = 0; i < order; i++)
            {
                A = lowerTriangle[i];
                B = upperTriangle[i];

                inner = Y[i];
                for (j = 0; j < i; j++)
                {
                    inner += A[j] * Y[j];
                }
                inner *= diagonal[i];

                result[i] = inner;
                for (j = 0; j < i; j++)
                {
                    result[j] += inner * B[j];
                }
            }
        }

        /// <summary>
        /// Solve a square Toeplitz system with a right-side matrix.
        /// </summary>
        /// <param name="Y">The right-side matrix</param>
        /// <returns>The solution matrix.</returns>
        /// <exception cref="ArgumentNullException">
        /// If <paramref name="Y"/> is <b>null</b>.
        /// </exception>
        /// <exception cref="NotConformableException">
        /// <para>The number of rows in <paramref name="Y"/> is not equal to the number of rows in the Toeplitz matrix.</para>
        /// </exception>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        /// <remarks>
        /// This member solves the linear system <b>TX</b> = <b>Y</b>, where <b>T</b> is
        /// a square Toeplitz matrix, <b>X</b> is the unknown solution matrix
        /// and <b>Y</b> is a known matrix.
        /// <para>
        /// The class implicitly decomposes the inverse Toeplitz matrix into a <b>UDL</b> factorisation
        /// using the Levinson algorithm, before calculating the solution vector.
        /// </para>
        /// </remarks>
        public Matrix Solve(Matrix Y)
        {
            // check parameters
            if (Y == null)
            {
                throw new System.ArgumentNullException("Y", Strings.NullParameterException);
            }

            if (order != Y.Rows)
            {
                throw new NotConformableException("Y", Strings.ParameterNotConformable);
            }

            if (IsSingular())
            {
                throw new SingularMatrixException();
            }
            Matrix result = MatrixBuilder.CreateMatrix(order, Y.Columns, MatrixType.Dense);
            Solve(Y, result);
            return result;
        }

        /// <summary>
        /// Solve a square Toeplitz system with a right-side matrix.
        /// </summary>
        /// <param name="Y">The right-side matrix</param>
        /// <param name="result">The solution matrix.</param>
        /// <exception cref="ArgumentNullException">
        /// If either <paramref name="Y"/> or <paramref name="result"/>is <b>null</b>.
        /// </exception>
        /// <exception cref="NotConformableException">
        /// <para>The number of rows in <paramref name="Y"/> is not equal to the number of rows in the Toeplitz matrix.</para>
        /// <para>- or -</para>
        /// <para>The number of rows in <paramref name="result"/> is not equal to the number of columns in the Toeplitz matrix.</para>
        /// <para>- or -</para>
        /// <para>The number of columns in <paramref name="result"/> is not equal to the number of columns in <paramref name="Y"/>.</para>
        /// </exception>
        /// <exception cref="SingularMatrixException">
        /// The Toeplitz matrix or one of the the leading sub-matrices is singular.
        /// </exception>
        /// <remarks>
        /// This member solves the linear system <b>TX</b> = <b>Y</b>, where <b>T</b> is
        /// a square Toeplitz matrix, <b>X</b> is the unknown solution matrix
        /// and <b>Y</b> is a known matrix.
        /// <para>
        /// The class implicitly decomposes the inverse Toeplitz matrix into a <b>UDL</b> factorisation
        /// using the Levinson algorithm, before calculating the solution vector.
        /// </para>
        /// </remarks>
        public void Solve(Matrix Y, Matrix result)
        {
            // check parameters
            if (Y == null)
            {
                throw new System.ArgumentNullException("Y", Strings.NullParameterException);
            }

            if (order != Y.Rows)
            {
                throw new NotConformableException("Y", Strings.ParameterNotConformable);
            }

            if (result == null)
            {
                throw new System.ArgumentNullException("result", Strings.NullParameterException);
            }

            if (order != result.Columns)
            {
                throw new NotConformableException("result", Strings.ParameterNotConformable);
            }

            if (Y.Columns != result.Columns)
            {
                throw new NotConformableException(Strings.ParametersNotConformable);
            }

            if (IsSingular())
            {
                throw new SingularMatrixException();
            }

            // allocate memory for solution
            double inner;
            double[] a, b;
            int i, j, l;
            double[] x = new double[order];
            double[] y = new double[order];

            for (l = 0; l < Y.Rows; l++)
            {
                // get right-side column
                for (i = 0; i < order; i++)
                {
                    y[i] = Y.ValueAt(i, l);
                }

                // solve left-side column
                for (i = 0; i < order; i++)
                {
                    a = lowerTriangle[i];
                    b = upperTriangle[i];

                    inner = y[i];
                    for (j = 0; j < i; j++)
                    {
                        inner += a[j] * y[j];
                    }
                    inner *= diagonal[i];

                    x[i] = inner;
                    for (j = 0; j < i; j++)
                    {
                        x[j] += inner * b[j];
                    }
                }

                // insert left-side column into the matrix
                result.SetColumn(l, x);
            }
        }

        #endregion Public Members

        #region Protected Members

        /// <summary>
        /// Delegates the actual computation to the correct Levinson subtype.
        /// </summary>
        protected override void InternalCompute()
        {
            int i, j, l;
            double Q, S;			// inner products
            double Ke, Kr;			// reflection coefficients
            double[] AP;			// reference to previous row in lower triangle
            double[] A;				// reference to current row in lower triangle
            double[] BP;			// reference to previous row in upper triangle
            double[] B;				// reference to current row in upper triangle

            // check if leading diagonal is zero
            if (leftColumn[0] == 0.0)
            {
                isSingular = true;
                return;
            }

            // setup zero order solution
            AP = lowerTriangle[0];
            AP[0] = 1.0;
            diagonal[0] = 1.0 / leftColumn[0];
            BP = upperTriangle[0];
            BP[0] = 1.0;

            for (i = 1; i < order; i++, AP = A, BP = B)
            {
                // calculate inner products
                Q = 0.0;
                for (j = 0, l = 1; j < i; j++, l++)
                {
                    Q += leftColumn[l] * AP[j];
                }

                S = 0.0;
                for (j = 0, l = 1; j < i; j++, l++)
                {
                    S += topRow[l] * BP[j];
                }

                // reflection coefficients
                Kr = -S * diagonal[i - 1];
                Ke = -Q * diagonal[i - 1];

                // update lower triangle
                A = lowerTriangle[i];
                A[0] = 0.0;

                Buffer.BlockCopy(AP, 0, A, Constants.SizeOfDouble, i * Constants.SizeOfDouble);
                for (j = 0, l = i - 1; j < i; j++, l--)
                {
                    A[j] += Ke * BP[l];
                }

                // update upper triangle
                B = upperTriangle[i];
                B[0] = 0.0;

                Buffer.BlockCopy(BP, 0, B, Constants.SizeOfDouble, i * Constants.SizeOfDouble);
                for (j = 0, l = i - 1; j < i; j++, l--)
                {
                    B[j] += Kr * AP[l];
                }

                // check for singular sub-matrix)
                if (Ke * Kr == 1.0)
                {
                    isSingular = true;
                    return;
                }

                // diagonal matrix
                diagonal[i] = diagonal[i - 1] / (1.0 - Ke * Kr);
            }
        }
        #endregion Protected Members

        #region Private Memebers
        private void CreateArrays()
        {
            // allocate memory for lower triangular matrix
            lowerTriangle = new double[order][];
            for (int i = 0; i < order; i++)
            {
                lowerTriangle[i] = new double[i + 1];
            }

            // allocate memory for diagonal
            diagonal = new double[order];

            // allocate memory for upper triangular matrix
            upperTriangle = new double[order][];
            for (int i = 0; i < order; i++)
            {
                upperTriangle[i] = new double[i + 1];
            }
        }
        #endregion Private Members
    }
}
